Ben Doran and Kyle Kimler

Compiled: August 13, 2021

ARBOL iteratively explores axes of variation in scRNAseq data by clustering and subclustering until variation between cells in subclusters becomes noise. The philosophy of ARBOL is that every axis of variation could be biologically meaningful, so each should be explored. Once every possibility is explored, curation and a statistical interrogation of the resolution are used to collapse clusters into the elemental transcriptomes of the dataset. ARBOL inherently builds a tree of subclustering events. As data is separated by the major axes of variation in each subset, further rounds capture less pronounced variables. This comes with some caveats: variation shared by all celltypes make up one of the major axes of variation in the first round of clustering. Celltypes can split up at the beginning, so the same splitting of e.g. B and T cells might happen further down in separate branches. Thus the resulting clustering tree is neither indicitave of true distances between end clusters nor a tree of unique groupings. We address this by manual curation of cell types, subtypes, and states, using literature-defined markers in the end clusters, an end-cluster naming method that captures the unique markers of the end clusters, and the calculation of a binary tree of these manually curated end clusters, providing a distance between the elemental transcriptomes of the dataset.

ARBOL was built to be a modular software. Currently it iterates through the default Seurat analysis pipeline, with two additions that enable automation:
- PCs used in clustering are chosen by 15 percent improvement of variance explained
- clustering parameters are chosen by maximization of silhouette size on a downsampled parameter scan

Normalization, dimensionality reduction, clustering, and plotting methods are set apart in the code in modularized functions to enable customization.

We saw a need for this software to address the fact that we were constantly manually subclustering scRNAseq data, and because other iterative clustering softwares like scPanoView and iterclust are tailor built for specific clustering methods outside of the Seurat paradigm.

This tutorial uses the following model datasets:
- pediatric FGID from the paper: 89849 cells x 22431 genes
- Nasal Polyp (Ordovas-Montanes et al 2018): 19196 cells x 22360 genes

FGID ARBOL

Reproducing figures for the 2021 paper

Running ARBOL

The github page describes ARBOL
ARBOL is called directly on a seurat object:

srobj <- readRDS("/path/to/full_seurat_object.rds")
endclustSrobjs <- GenTieredclusters(srobj,
                           saveSROBJdir = "~/tieredoutput/srobjs",
                           figdir = "~/tieredoutput/figdir",
                           SaveEndNamesDir = "~/tieredoutput/endclusts")

Specifying output directories saves QC statistics



clustering decision metrics



differential expression results among clusters (and visualizations) at each level in ‘figdir’

##            p_val avg_logFC pct.1 pct.2     p_val_adj cluster  gene
## 1: 7.468466e-307  2.148874 0.622 0.061 1.695342e-302       0  CD3D
## 2: 1.898258e-290  2.127852 0.664 0.128 4.309047e-286       0  TRAC
## 3: 3.575767e-214  1.774795 0.479 0.045 8.116990e-210       0  CD3E
## 4: 5.584884e-210  1.798812 0.644 0.225 1.267769e-205       0  IL32
## 5: 9.468750e-209  2.275919 0.472 0.055 2.149406e-204       0 KLRB1
## 6: 1.609778e-208  1.882106 0.480 0.044 3.654196e-204       0   CD2


and UMAP visualizations at each tier.



Finally, GenTieredClusters() also saves subset seurat objects for each cluster in ‘srobjs’, and saves files with lists of cells per end-cluster in ‘endclusts’.

Iterative clustering is enabled by automation of two key parameter choices. Principle components are chosen using two heuristics: when more than 500 cells are present, PCs are included if they explain 15% more variance than the following component. When less than 500 cells, Seurat’s Jackstraw method is used to calculate significant principle components. Clustering is performed automatically by performing a parameter scan on a downsampled dataset using Seurat’s built-in Louvain clustering. Silhouette measures the ratio of intra-cluster distance to inter-cluster distance, where a high score means highly distinct clusters. For stages where we were clustering more than 500 cells, a randomized subsample of N cells / 10 was used in the parameter scane.



Working with ARBOL results

I prefer to load tiered cluster data into a dataframe, which I can use directly for comparisons or which I can reattach to the full seurat object’s metadata for annotation, dendrogram production, and other analysis.

LoadTiersAsDF <- function(folder='./tieredoutput/endclusts',maxtiers=10) {
  tierfiles <- list.files(folder,full.names=TRUE,recursive=TRUE)
  sample_strings <- sub('\\..*$','',basename(tierfiles))
  sample_strings <- sub('_T0C0_', '',sample_strings)
  tiers <- map2(tierfiles, sample_strings, ~fread(.x,header=FALSE) %>% mutate(id = .y))
  tiers <- rbindlist(tiers) %>% data.frame
  tiers$id <- gsub('_','.',tiers$id)
  tiers <- tiers %>% separate(id,into=paste0('tier',1:maxtiers),sep='\\.',remove=FALSE)
  tiers <- tiers %>% rename(CellID=V1)
  tiers <- tiers %>% data.frame %>% rename(tierNident=id)
  return(tiers)
}
fgtiers <- LoadTiersAsDF('~/PREDICT_3p_Paper/FGID/tieredoutput/endclusts')
print(head(fgtiers))
##                                   CellID               tierNident tier1 tier2
## 1 p049_T0D_ILE_LPS_3p_AAGACCTTCCCATTTA-1 T0C0.T1C4.T2C0.T3C5.T4C3  T0C0  T1C4
## 2 p049_T0D_ILE_LPS_3p_ACGGGCTTCAAACCAC-1 T0C0.T1C4.T2C0.T3C5.T4C3  T0C0  T1C4
## 3 p049_T0D_ILE_LPS_3p_ACTGAGTAGTTCGATC-1 T0C0.T1C4.T2C0.T3C5.T4C3  T0C0  T1C4
## 4 p049_T0D_ILE_LPS_3p_AGCTTGACAGCCTGTG-1 T0C0.T1C4.T2C0.T3C5.T4C3  T0C0  T1C4
## 5 p049_T0D_ILE_LPS_3p_AGTTGGTTCACTGGGC-1 T0C0.T1C4.T2C0.T3C5.T4C3  T0C0  T1C4
## 6 p049_T0D_ILE_LPS_3p_CATGGCGAGACGCTTT-1 T0C0.T1C4.T2C0.T3C5.T4C3  T0C0  T1C4
##   tier3 tier4 tier5 tier6 tier7 tier8 tier9 tier10
## 1  T2C0  T3C5  T4C3  <NA>  <NA>  <NA>  <NA>   <NA>
## 2  T2C0  T3C5  T4C3  <NA>  <NA>  <NA>  <NA>   <NA>
## 3  T2C0  T3C5  T4C3  <NA>  <NA>  <NA>  <NA>   <NA>
## 4  T2C0  T3C5  T4C3  <NA>  <NA>  <NA>  <NA>   <NA>
## 5  T2C0  T3C5  T4C3  <NA>  <NA>  <NA>  <NA>   <NA>
## 6  T2C0  T3C5  T4C3  <NA>  <NA>  <NA>  <NA>   <NA>

Add ARBOL annotations to the seurat object

srobj <- readRDS('~/PREDICT_3p_Paper/FGID/srobj.rds')
srobj@meta.data <- srobj@meta.data %>% cbind(as.data.frame(str_split_fixed(srobj@meta.data$CellID,"_",n=6))) %>% dplyr::rename(patient = V1, time_point = V2, tissue = V3, sample_site = V4, seqDirection=V5,barcode=V6)

srobj@meta.data <- srobj@meta.data %>% left_join(fgtiers,by="CellID")

print(head(srobj@meta.data))
##                                   CellID          curatedname patient
## 1 p009_T0D_ILE_LPS_3p_AAACCTGAGACTGTAA-1   T/NK/ILC.MAF.RPS26    p009
## 2 p009_T0D_ILE_LPS_3p_AAACCTGAGTAGCCGA-1   T/NK/ILC.CCR7.SELL    p009
## 3 p009_T0D_ILE_LPS_3p_AAACCTGAGTAGGCCA-1 Endth/Cap.PLVAP.FLT1    p009
## 4 p009_T0D_ILE_LPS_3p_AAACCTGAGTCAAGGC-1   T/NK/ILC.CCR7.SELL    p009
## 5 p009_T0D_ILE_LPS_3p_AAACCTGCAACGATGG-1         B.IGHD.FCER2    p009
## 6 p009_T0D_ILE_LPS_3p_AAACCTGCAATCCGAT-1   T/NK/ILC.CCR7.SELL    p009
##   time_point tissue sample_site seqDirection            barcode
## 1        T0D    ILE         LPS           3p AAACCTGAGACTGTAA-1
## 2        T0D    ILE         LPS           3p AAACCTGAGTAGCCGA-1
## 3        T0D    ILE         LPS           3p AAACCTGAGTAGGCCA-1
## 4        T0D    ILE         LPS           3p AAACCTGAGTCAAGGC-1
## 5        T0D    ILE         LPS           3p AAACCTGCAACGATGG-1
## 6        T0D    ILE         LPS           3p AAACCTGCAATCCGAT-1
##                 tierNident tier1 tier2 tier3 tier4 tier5 tier6 tier7 tier8
## 1      T0C0.T1C0.T2C1.T3C0  T0C0  T1C0  T2C1  T3C0  <NA>  <NA>  <NA>  <NA>
## 2      T0C0.T1C0.T2C1.T3C1  T0C0  T1C0  T2C1  T3C1  <NA>  <NA>  <NA>  <NA>
## 3           T0C0.T1C8.T2C1  T0C0  T1C8  T2C1  <NA>  <NA>  <NA>  <NA>  <NA>
## 4      T0C0.T1C0.T2C1.T3C1  T0C0  T1C0  T2C1  T3C1  <NA>  <NA>  <NA>  <NA>
## 5 T0C0.T1C1.T2C0.T3C1.T4C3  T0C0  T1C1  T2C0  T3C1  T4C3  <NA>  <NA>  <NA>
## 6      T0C0.T1C0.T2C1.T3C1  T0C0  T1C0  T2C1  T3C1  <NA>  <NA>  <NA>  <NA>
##   tier9 tier10
## 1  <NA>   <NA>
## 2  <NA>   <NA>
## 3  <NA>   <NA>
## 4  <NA>   <NA>
## 5  <NA>   <NA>
## 6  <NA>   <NA>

End cluster diversity analysis

In the paper, we used a diversity metric, the Gini-Simpson’s Index, to measure diversity of patients in our end clusters to determine if the cluster was powered for differential analysis. Gini-Simpson’s index is calculated as

\(diversity = 1 - \lambda = \sum_{i=1}^{R} p_{i}^{2}\)

R is the total number of patients in the dataset while \(p_{i}\) is the proportional abundance of each patient in the end clusters. lambda is subtracted from 1 because in datasets with fewer patients, lambda will be inflated.

1 indicates equal representation of all patients within a subset and 0 indicates a completely patient specific subset

tierNdiversity <- srobj@meta.data %>% group_by(tierNident) %>% 
        summarize(diversity = 1 - sum(((table(patient)) / (length(patient)))^2)) 
curatedNamediversity <- srobj@meta.data %>% group_by(curatedname) %>% 
        summarize(diversity = 1 - sum(((table(patient)) / (length(patient)))^2))

We used patient diversity among end clusters to inform manual curation

s1 <- ggplot(tierNdiversity,aes(x=diversity))+ geom_histogram() + ggtitle("raw ARBOL results")
s2 <- ggplot(curatedNamediversity,aes(x=diversity))+ geom_histogram() + ggtitle("after manual curation")
s1+s2

Building a tree

meta <- srobj@meta.data
meta$pathString <- meta$tierNident
meta$pathString <- gsub('\\.','/',meta$pathString)
ARBOLtree <- as.Node(meta) 
print(ARBOLtree,'majortype','subtype','groundtruth','curatedname',limit=20)
##                                   levelName majortype subtype groundtruth
## 1  T0C0                                            NA      NA          NA
## 2   ¦--T1C0                                        NA      NA          NA
## 3   ¦   ¦--T2C1                                    NA      NA          NA
## 4   ¦   ¦   ¦--T3C0                                NA      NA          NA
## 5   ¦   ¦   °--T3C1                                NA      NA          NA
## 6   ¦   ¦--T2C0                                    NA      NA          NA
## 7   ¦   ¦   ¦--T3C0                                NA      NA          NA
## 8   ¦   ¦   ¦   ¦--T4C0                            NA      NA          NA
## 9   ¦   ¦   ¦   ¦   ¦--T5C8                        NA      NA          NA
## 10  ¦   ¦   ¦   ¦   ¦--T5C7                        NA      NA          NA
## 11  ¦   ¦   ¦   ¦   ¦--T5C0                        NA      NA          NA
## 12  ¦   ¦   ¦   ¦   ¦--T5C5                        NA      NA          NA
## 13  ¦   ¦   ¦   ¦   ¦--T5C6                        NA      NA          NA
## 14  ¦   ¦   ¦   ¦   ¦--T5C1                        NA      NA          NA
## 15  ¦   ¦   ¦   ¦   ¦--T5C9                        NA      NA          NA
## 16  ¦   ¦   ¦   ¦   ¦--T5C3                        NA      NA          NA
## 17  ¦   ¦   ¦   ¦   ¦   ¦--T6C0                    NA      NA          NA
## 18  ¦   ¦   ¦   ¦   ¦   ¦--T6C1                    NA      NA          NA
## 19  ¦   ¦   ¦   ¦   ¦   °--T6C2                    NA      NA          NA
## 20  ¦   ¦   ¦   ¦   °--... 2 nodes w/ 0 sub        NA      NA          NA
## 21  ¦   ¦   ¦   °--... 3 nodes w/ 5 sub            NA      NA          NA
## 22  ¦   ¦   °--... 3 nodes w/ 23 sub               NA      NA          NA
## 23  ¦   °--... 1 nodes w/ 26 sub                   NA      NA          NA
## 24  °--... 11 nodes w/ 403 sub                     NA      NA          NA
##           curatedname
## 1                    
## 2                    
## 3                    
## 4  T/NK/ILC.MAF.RPS26
## 5  T/NK/ILC.CCR7.SELL
## 6                    
## 7                    
## 8                    
## 9         T.RORA.CCR6
## 10        T.RORA.CCR6
## 11        T.RORA.CCR6
## 12        T.RORA.CCR6
## 13        T.RORA.CCR6
## 14        T.RORA.CCR6
## 15        T.RORA.CCR6
## 16                   
## 17        T.RORA.CCR6
## 18        T.RORA.CCR6
## 19        T.RORA.CCR6
## 20                   
## 21                   
## 22                   
## 23                   
## 24
ARBOLphylo <- as.phylo(ARBOLtree)

#For curatednames
#ARBOLphylo$tip.label <- ARBOLtree$Get('curatedname')[!is.na(ARBOLtree$Get('curatedname'))]
#For tierNident
ARBOLphylo$tip.label <- ARBOLtree$Get('tierNident')[!is.na(ARBOLtree$Get('tierNident'))]
ggtree(ARBOLphylo,branch.length='edge.length',layout='circular')

Display Gini-Simpson’s index next to the tree, in the paper’s style

plotdendrogram <- function(phyloObj,title) {
    dd.row <- as.dendrogram(phyloObj)
    ddata_x <- dendro_data(dd.row)
    ddata_L <- label(ddata_x)
    ddata_L$label <- str_replace_all(ddata_L$label, "−", "-")
    #For majortype
    #ddata_L$group <- gsub("\\..*$", "", ddata_L$label)
    #For subtype
    ddata_L$group <- gsub("^([^.]+)\\.([^.]+)\\.(.+)$","\\2",ddata_L$label)
    ddata_n <- segment(ddata_x) %>% filter(yend!=0 & yend!=y)
    celltype.colors <- primary.colors(length(unique(ddata_L$group)))
    names(celltype.colors) <- unique(ddata_L$group)

    ddata_L$count = as.numeric(table(factor(srobj$tierNident))[ARBOLphylo$tip.label])

    ddata_L$diversity = srobj@meta.data %>% group_by(tierNident) %>% 
        summarize(diversity = 1 - sum(((table(patient)) / (length(patient)))^2)) %>%
        column_to_rownames('tierNident') %>%
        .[phyloObj$tip.label, ]

    y.scl = .05
    plt <- ggplot(segment(ddata_x)) +
        geom_segment(aes(x=x, y=y*y.scl, xend=xend, yend=yend*y.scl)) + 
        geom_text(data=ddata_L, angle = 0, hjust = 0, size = 1.5,
                  aes(label=label, x=x, y=-.4, color=group)) + 
        geom_text(data=ddata_L, angle = 0, hjust = 1, size = 1.5,
                  aes(label=count, x=x, y=-0.39),color="black") +
        geom_tile(data=ddata_L, aes(x = x, y=-0.11, fill=diversity), height=0.15) +
        scale_color_manual(values=celltype.colors) +
        scale_fill_gradient(low = "white", high = "black", limits=c(0,1)) +
        coord_flip(clip = 'off') + 
        scale_y_reverse(limits = c(max(segment(ddata_x)$y*y.scl),-4)) +
                           #breaks = c(1, 0.75, 0.5, 0.25, 0),
                           #labels = c(1, 0.75, 0.5, 0.25, 0)) +
        scale_x_reverse() +
        theme_classic() + 
        theme(legend.position = "bottom",
              axis.text.x = element_blank(),
              axis.ticks.x = element_blank()) +
        ylab("Split Distance") +
        xlab("Clusters") +
        ggtitle(title)


    return(plt)
}
p1 <- plotdendrogram(ARBOLphylo, title='raw ARBOL')
p1

Annotation

In the paper, tiers 1 and 2 were well annotated by comparing Wilcoxon 1 v all results to markers in the literature. End clusters were named based on the tier1.tier2 combination they belonged to, followed by 2 genes either manually chosen or chosen on a heuristic rank score among its tier2 group. The rank score is calculated as follows on the results of a 1 vs. all wilcoxon rank sum test:

\[-\log (p.adj+(1*10^{-310})) * \log{foldchange} * \small\cfrac{pct.1}{pct.2+1*10^{-300}}\]

We also created a shiny app for understanding how the weighting the components of the rank score affects gene choice

We can view the raw ARBOL with curated names to get an idea of how we collapsed clusters

curatedNamediversity <- srobj@meta.data %>% group_by(curatedname) %>% 
        summarize(diversity = 1 - sum(((table(patient)) / (length(patient)))^2)) 
ARBOLphylo$tip.label <- ARBOLtree$Get('curatedname')[!is.na(ARBOLtree$Get('curatedname'))]
p2 <- plotdendrogram(as.phylo(ARBOLphylo),"raw ARBOL with curated names")
p2

It is important to note that ARBOL can split major cell types in the earliest stages, so that the final tree may misrepresent the relationships of the end clusters. To reconstruct those relationships, we hierarchically clustered euclidean distances between gene medians among end clusters, building a binary tree

Binary tree creation

While curating the end clusters, we wanted to better understand the distance between clusters. In the atlas figures 3 and 4, we performed hierarchical clustering of curated end clusters (not to be confused with the raw ARBOL: hclus produces a binary tree). Hierarchical clustering was performed on median expression of a subset of highly DE genes (approximately 4500) using complete linkage and distance calculated with Pearson correlation. DE was performed by pairwise expression tests of end clusters.

Idents(srobj) <- srobj@meta.data$curatedname
#markers are selected by pairwise expression tests
de.df.sct <- lapply(Idents(srobj), function(x) {return(FindMarkers(srobj,ident.1=x))})
de.df.sct <- bind_rows(de.df.sct)
genes.use <- de.df.sct$gene

scaled.data = Matrix::Matrix(srobj@assays$SCT@counts[de.df.sct$gene,])
scaled.data.t = t(scaled.data)
scaled.data.t = as.data.frame(scaled.data.t)
scaled.data.t$cellsubsets = srobj$curatedname

agg.clst.cntrs = aggregate(scaled.data.t[, -which(colnames(scaled.data.t) == "cellsubsets")],
          list(scaled.data.t$cellsubsets), median)

colnames(agg.clst.cntrs) <- c("cellsubsets", colnames(agg.clst.cntrs)[2:length(colnames(agg.clst.cntrs))])

clst.cntrs <- agg.clst.cntrs %>% 
    column_to_rownames("cellsubsets") %>%
    t %>% as.data.frame

The package pvclust was used for hierarchical clustering, bootstrapping sampling genes. Here I use one bootstrap due to computational limitations

library(pvclust)
result <- pvclust(clst.cntrs, method.dist="cor", method.hclust="complete", nboot=1, parallel=TRUE)
## Creating a temporary cluster...done:
## socket cluster with 11 nodes on host 'localhost'
## Warning in pvclust.parallel(cl = cl, data = data, method.hclust =
## method.hclust, : nboot is too small for cluster size. non-parallel version is
## executed
## Bootstrap (r = 0.5)... Done.
## Bootstrap (r = 0.6)... Done.
## Bootstrap (r = 0.7)... Done.
## Bootstrap (r = 0.8)... Done.
## Bootstrap (r = 0.9)... Done.
## Bootstrap (r = 1.0)... Done.
## Bootstrap (r = 1.1)... Done.
## Bootstrap (r = 1.2)... Done.
## Bootstrap (r = 1.3)... Done.
## Bootstrap (r = 1.4)... Done.
p3

Nasal Polyp ARBOL

Here is another example of an ARBOL built on another scRNAseq dataset of nasal polyps.

## Warning: Expected 10 pieces. Missing pieces filled with `NA` in 19196 rows [1,
## 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, ...].
tiers <- LoadTiersAsDF(folder = '~/tieredoutput/endclusts')
nasalp1

nasalp2

nasalp3

## [1] "Thanks for reading!"
sessionInfo()
## R version 4.1.0 (2021-05-18)
## Platform: x86_64-apple-darwin17.0 (64-bit)
## Running under: macOS Big Sur 10.16
## 
## Matrix products: default
## BLAS:   /Library/Frameworks/R.framework/Versions/4.1/Resources/lib/libRblas.dylib
## LAPACK: /Library/Frameworks/R.framework/Versions/4.1/Resources/lib/libRlapack.dylib
## 
## locale:
## [1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8
## 
## attached base packages:
## [1] grid      stats     graphics  grDevices utils     datasets  methods  
## [8] base     
## 
## other attached packages:
##  [1] phylogram_2.1.0    dendextend_1.15.1  cowplot_1.1.1      phyloseq_1.36.0   
##  [5] rdist_0.0.5        ape_5.5            ggtree_3.0.2       treeio_1.16.1     
##  [9] data.tree_1.0.0    Matrix_1.3-4       ggpubr_0.4.0       glmGamPoi_1.4.0   
## [13] vegan_2.5-7        lattice_0.20-44    permute_0.9-5      pheatmap_1.0.12   
## [17] ggdendro_0.1.22    pvclust_2.2-0      colorRamps_2.3     networkD3_0.4     
## [21] viridis_0.6.1      viridisLite_0.4.0  SeuratObject_4.0.2 Seurat_4.0.1      
## [25] data.table_1.14.0  forcats_0.5.1      stringr_1.4.0      dplyr_1.0.7       
## [29] purrr_0.3.4        readr_2.0.1        tidyr_1.1.3        tibble_3.1.3      
## [33] ggplot2_3.3.5      tidyverse_1.3.1    rmarkdown_2.10    
## 
## loaded via a namespace (and not attached):
##   [1] utf8_1.2.2                  reticulate_1.20            
##   [3] tidyselect_1.1.1            htmlwidgets_1.5.3          
##   [5] Rtsne_0.15                  munsell_0.5.0              
##   [7] codetools_0.2-18            ica_1.0-2                  
##   [9] future_1.21.0               miniUI_0.1.1.1             
##  [11] withr_2.4.2                 colorspace_2.0-2           
##  [13] Biobase_2.52.0              highr_0.9                  
##  [15] knitr_1.33                  rstudioapi_0.13            
##  [17] stats4_4.1.0                ROCR_1.0-11                
##  [19] ggsignif_0.6.2              tensor_1.5                 
##  [21] listenv_0.8.0               MatrixGenerics_1.4.2       
##  [23] labeling_0.4.2              GenomeInfoDbData_1.2.6     
##  [25] polyclip_1.10-0             bit64_4.0.5                
##  [27] farver_2.1.0                rhdf5_2.36.0               
##  [29] parallelly_1.27.0           vctrs_0.3.8                
##  [31] generics_0.1.0              xfun_0.25                  
##  [33] R6_2.5.0                    GenomeInfoDb_1.28.1        
##  [35] rhdf5filters_1.4.0          bitops_1.0-7               
##  [37] spatstat.utils_2.2-0        DelayedArray_0.18.0        
##  [39] assertthat_0.2.1            vroom_1.5.4                
##  [41] promises_1.2.0.1            scales_1.1.1               
##  [43] gtable_0.3.0                globals_0.14.0             
##  [45] goftest_1.2-2               rlang_0.4.11               
##  [47] splines_4.1.0               rstatix_0.7.0              
##  [49] lazyeval_0.2.2              spatstat.geom_2.2-2        
##  [51] broom_0.7.9                 BiocManager_1.30.16        
##  [53] yaml_2.2.1                  reshape2_1.4.4             
##  [55] abind_1.4-5                 modelr_0.1.8               
##  [57] backports_1.2.1             httpuv_1.6.1               
##  [59] tools_4.1.0                 ellipsis_0.3.2             
##  [61] spatstat.core_2.3-0         biomformat_1.20.0          
##  [63] jquerylib_0.1.4             RColorBrewer_1.1-2         
##  [65] BiocGenerics_0.38.0         ggridges_0.5.3             
##  [67] Rcpp_1.0.7                  plyr_1.8.6                 
##  [69] zlibbioc_1.38.0             RCurl_1.98-1.3             
##  [71] rpart_4.1-15                deldir_0.2-10              
##  [73] pbapply_1.4-3               S4Vectors_0.30.0           
##  [75] zoo_1.8-9                   SummarizedExperiment_1.22.0
##  [77] haven_2.4.3                 ggrepel_0.9.1              
##  [79] cluster_2.1.2               fs_1.5.0                   
##  [81] magrittr_2.0.1              scattermore_0.7            
##  [83] openxlsx_4.2.4              lmtest_0.9-38              
##  [85] reprex_2.0.1                RANN_2.6.1                 
##  [87] fitdistrplus_1.1-5          matrixStats_0.60.0         
##  [89] hms_1.1.0                   patchwork_1.1.1            
##  [91] mime_0.11                   evaluate_0.14              
##  [93] xtable_1.8-4                rio_0.5.27                 
##  [95] readxl_1.3.1                IRanges_2.26.0             
##  [97] gridExtra_2.3               compiler_4.1.0             
##  [99] KernSmooth_2.23-20          crayon_1.4.1               
## [101] htmltools_0.5.1.1           mgcv_1.8-36                
## [103] later_1.2.0                 tzdb_0.1.2                 
## [105] aplot_0.0.6                 lubridate_1.7.10           
## [107] DBI_1.1.1                   dbplyr_2.1.1               
## [109] MASS_7.3-54                 ade4_1.7-17                
## [111] car_3.0-11                  cli_3.0.1                  
## [113] parallel_4.1.0              igraph_1.2.6               
## [115] GenomicRanges_1.44.0        pkgconfig_2.0.3            
## [117] rvcheck_0.1.8               foreign_0.8-81             
## [119] plotly_4.9.4.1              spatstat.sparse_2.0-0      
## [121] foreach_1.5.1               xml2_1.3.2                 
## [123] bslib_0.2.5.1               multtest_2.48.0            
## [125] XVector_0.32.0              rvest_1.0.1                
## [127] digest_0.6.27               sctransform_0.3.2          
## [129] RcppAnnoy_0.0.19            Biostrings_2.60.2          
## [131] spatstat.data_2.1-0         cellranger_1.1.0           
## [133] leiden_0.3.9                tidytree_0.3.4             
## [135] uwot_0.1.10                 curl_4.3.2                 
## [137] shiny_1.6.0                 lifecycle_1.0.0            
## [139] nlme_3.1-152                jsonlite_1.7.2             
## [141] Rhdf5lib_1.14.2             carData_3.0-4              
## [143] fansi_0.5.0                 pillar_1.6.2               
## [145] fastmap_1.1.0               httr_1.4.2                 
## [147] survival_3.2-11             glue_1.4.2                 
## [149] zip_2.2.0                   iterators_1.0.13           
## [151] png_0.1-7                   bit_4.0.4                  
## [153] stringi_1.7.3               sass_0.4.0                 
## [155] irlba_2.3.3                 future.apply_1.8.1